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ABSTRACT 

Recent observations of gaps and non-axisymmetric features in the dust distributions of transition 
disks have been interpreted as evidence of embedded massive protoplanets. However, comparing 
the predictions of planet-disk interaction models to the observed features has shown far from perfect 
agreement. This may be due to the strong approximations used for the predictions. For example, spi¬ 
ral arm fitting typically uses results that are based on low-mass planets in an isothermal gas. In this 
work, we describe two-dimensional, global, hydrodynamical simulations of disks with embedded 
protoplanets, with and without the assumption of local isothermality, for a range of planet-to-star 
mass ratios 1-10 Mj for a 1 M© star. We use the PENCIL CODE in polar coordinates for our models. 
We find that the inner and outer spiral wakes of massive protoplanets (M > 5 Mj) produce significant 
shock heating that can trigger buoyant instabilities. These drive sustained turbulence throughout the 
disk when they occur. The strength of this effect depends strongly on the mass of the planet and the 
thermal relaxation timescale; for a 10 Mj planet embedded in a thin, purely adiabatic disk, the spi¬ 
rals, gaps, and vortices typically associated with planet-disk interactions are disrupted. We find that 
the effect is only weakly dependent on the initial radial temperature profile. The spirals that form 
in disks heated by the effects we have described may fit the spiral structures observed in transition 
disks better than the spirals predicted by linear isothermal theory. 

Subject headings: hydrodynamics — planet-disk interactions — planets and satellites: formation — 
protoplanetary disks — shock waves — turbulence 


1. INTRODUCTION 

Decades of analytical and numerical work (Pa- 
paloizou & Lin 1984; Lin & Papaloizou 1986a,b; Nelson 
et al. 2000; Masset & Snellgrove 2001; Paardekooper & 
Mellema 2004; Quillen et al. 2004; de Val-Borro et al. 
2006; Klahr & Kley 2006; Lyra et al. 2009; Zhu et al. 
2011; Kley et al. 2012; Kley & Nelson 2012) have estab¬ 
lished that massive (> IMj) protoplanets generate ob¬ 
servable gaps and other features in their host disks that 
can in principle be resolved by ALMA (see, e.g., the 
review of Wolf et al. 2012). These models show that 
the gravitationally-induced tides of protoplanets with 
masses exceeding Saturn's clear axisymmetric gaps in 
the disk. These generate long-lived, non-axisymmetric 
vortices due to the Rossby wave instability (RWI) at the 
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gap edges (Lovelace & Hohlfeld 1978; Toomre 1981; Pa¬ 
paloizou & Pringle 1984, 1985; Hawley 1987; Lovelace 
et al. 1999). Non-axisymmetric dust clouds have been 
observed in transition disks at sub-mm wavelengths 
(Oppenheimer et al. 2008; Brown et al. 2009; Casassus 
et al. 2012; Isella et al. 2013; Casassus et al. 2013; van der 
Marel et al. 2013), features that are usually interpreted 
as vortices induced by embedded protoplanets. Spi¬ 
ral structures, one of the hallmarks of planet-disk inter¬ 
actions (Goldreich & Tremaine 1979, Ogilvie & Lubow 
2002, Rafikov 2002), are also seen in several transition 
disks in polarized scattered light (Muto et al. 2012, 
Garufi et al. 2013, Avenhaus et al. 2014, Currie et al. 
2014). 

However, because the putative corresponding plan¬ 
ets are not directly observed, sufficiently detailed dy¬ 
namical and radiative transfer models must be consis¬ 
tent with observations, and all other possible ways to 
produce the observed features must be ruled out. For 
example, there are several mechanisms that are known 
to generate vortices without a planet in the disk. The 
RWI can occur at dead zone boundaries (Varniere & Tag¬ 
ger 2006; Lyra et al. 2008), whether at the inner edge of 
the dead zone (Lyra & Mac Low 2012), or at the outer 
one (Lyra et al. 2015), despite the transition in resistiv¬ 
ity in the outer disk being smooth (Dzyurkevich et al. 
2013). In regions without the steep vortensity gradi¬ 
ent required for the RWI, a vortex could be maintained 
by baroclinic feedback as long as (1) a non-zero entropy 
gradient exists (Klahr & Bodenheimer 2003; Klahr 2004), 
(2) the thermal time is finite (Petersen et al. 2007a,b; 
Lesur & Papaloizou 2010; Raettig et al. 2013), and (3) 
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magnetization is absent (Lyra & Klahr 2011). Such con¬ 
ditions may or may not occur in the outer regions of pro¬ 
toplanetary disks. 

There are also difficulties in interpreting gaps as un¬ 
ambiguous signposts of planet-disk interactions. The 
transition disk surrounding HD 100546, for instance, 
contains a gap, but the gap edge is far too smooth to 
have been caused by a planet (Mulders et al. 2013). In 
optically thin disks, Lyra & Kuchner (2013) find that the 
combination of photoelectric heating and dust trapping 
can lead to the production of sharp rings similar to those 
commonly attributed to gravitational shepherding by 
planets. Furthermore, the magnetorotational instabil¬ 
ity (MRI; Balbus & Hawley 1991) could lead to gaps in 
the dust distribution by inducing zonal flows that are 
triggered by magnetic flux accumulation in turbulence 
(Lyra et al. 2008; Johansen et al. 2009; Simon et al. 2012; 
Kunz & Lesur 2013; Flock et al. 2015). 

Finally, a one-to-one correlation between planets and 
spirals is also absent. Shear causes any density wave 
in a differentially rotating disk to propagate as a spiral. 
Therefore, the observed spiral features could easily be 
the result of disk turbulence. For example, Lyra et al. 
(2015) show that waves from the MRI-active zone prop¬ 
agate into the MRI-dead zone as a coherent spiral. In 
fact, it has been difficult to explain the spiral features ob¬ 
served in transition disk as unequivocally due to plan¬ 
ets. The observed spirals have, in general, very large 
pitch angles, which is inconsistent with the background 
disk temperature if linear spiral density wave theory for 
an isothermal gas is assumed (Rafikov 2002; Muto et al. 
2012). Currie et al. (2014) present a remarkable spiral 
feature in the disk around HD 100546, with little polar¬ 
ization, that can only be fit with an effective disk tem¬ 
perature of ^1000 K, substantially higher than the ex¬ 
pected gas temperatures. Benisty et al. (2015) show that 
fitting a model to the spiral in the disk around MWC 758 
requires a large disk aspect ratio around the spiral fea¬ 
tures, corresponding to 300 K at 55AU, whereas the 
background gas is at much colder temperatures, 50 K. 
Juhasz et al. (2014) show that a local increase in pres¬ 
sure scale height by at least 20% would be required to 
reproduce observations of multiple disks in polarized 
scattered light. 

The high temperatures implied by these spiral fea¬ 
tures motivate us to abandon the idea of local isother- 
mality, and instead to entertain potential effects due 
to the departure from barotropic conditions. Although 
numerical simulations of planet-disk interactions have 
been performed for decades, a search of the literature re¬ 
veals that most previous works have either made use of 
the locally isothermal approximation or examined plan¬ 
etary masses of no more than ~l-5Mj. As some of 
the candidate planets in transition disks are very mas¬ 
sive (> 10 Mj), we also examine that region of the pa¬ 
rameter space. We find that the wake generated by the 
most massive planets has relative Mach numbers above 
unity compared to the quiescent disk. The resulting 
shocks heat the surrounding gas, effectively converting 
the planet's gravitational potential energy into a power¬ 
ful extra energy source for disk heating. For cases with 
sufficiently long cooling times (and in which the thin- 
disk approximation holds), a radial buoyant instability 


occurs, leading to turbulence around the planet and dis¬ 
ruption of the usual wake features, i.e., spirals, gaps, 
and vortices. 

This paper is the first in a series in which we com¬ 
pare the behavior of planets embedded in disks that 
are evolved with and without the assumption of local 
isothermality. The next paper in the series will extend 
the two-dimensional (2D) study presented here to in¬ 
clude three-dimensional disks. In Sect. 2, we present 
our numerical model model. We discuss the simulation 
results in Sect. 3, and conclude in Sect. 4. 

2. THE MODEL 
2.1. Governing equations 

In this work, we perform two-dimensional global hy- 
drodynamical simulations of gas disks with an embed¬ 
ded massive planet to explore the effects of shocks that 
result from planet-disk interactions in non-barotropic 
disks. We primarily use the PENCIL CODE^, a high- 
order finite-difference grid hydrodynamics code. An in¬ 
dependent code, BoxzyHydro, is also used to validate 
the main results, as described in Sect. 3.3. For the Pencil 
simulations, the gas density is evolved using the conti¬ 
nuity equation 


Dr 


= -r V • u 


and the equation of motion 


( 1 ) 


^ =-l;VP-V#+1;V (2) 

Dt E E ^ ^ ' 

where 27, u, and P are the surface density, velocity, 
and pressure of the gas, respectively, and ^ represents 
the gravitational potential contributions of the star and 
planet (we have ignored the effects of gas self-gravity). 
The operator 


D 

Dt 


dt 


+ !/• V 


(3) 


represents the advective derivative. The last term in the 
equation of motion is the acceleration due to shock vis¬ 
cosity. Assuming that the shock viscosity takes the 
form of a bulk viscosity, the stress tensor is 


C// = i^sh-^^!7 V-M, (4) 

The shock viscosity must be a localized function, en¬ 
suring that sufficient energy is dissipated in regions of 
the flow where shocks occur to satisfy the Rankine- 
Hugoniot jump conditions, while more quiescent re¬ 
gions are left unaffected. We take 


^Sh = Csh ^niax[(-V • r/)+]y [min(Ar,rA^)]^. (5) 

The viscosity is proportional to the maximum of posi¬ 
tive flow convergence, as evaluated over three grid cells 
in each direction for a total of 9 zones in 2D (the given 
cell plus its immediate neighbors). The angled brackets 


^ See http://pencil-code.googlecode.com 
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represent a quadratic smoothing function that smooths 
the divergence over seven zones in each direction, with 
weights (l,6,15,20,15,6,l)/64. The result is then scaled 
by the square of the smallest grid spacing. In the above 
equation, r is the cylindrical radius and (p the azimuthal 
angle; Ar and Acp give the numerical grid spacing. The 
quantity Cgh is a constant defining the strength of the 
shock viscosity, set to unity in our simulations. Varying 
this coefficient changes the area over which the added 
entropy is distributed but not the total amount of extra 
energy, which still must satisfy the jump conditions. We 
refer to Cg^ as the shock viscosity coefficient (for more 
details see Haugen et al. 2004). 

To follow the thermodynamic evolution of the disk, 
we assume the ideal gas law as the equation of state 

P = Scl/j, ( 6 ) 

where Cs is the adiabatic sound speed and 7 is the adi¬ 
abatic index. The temperature and sound speed are re¬ 
lated according to 


T = c2/[cp(7-l)], (7) 

where Cp = jCy is the specific heat at constant pressure, 
and Cy the specific heat at constant volume. They are 
related to the universal gas constant R by 

R/]i = Cp-Cy, 

where ]i is the mean molecular weight of the gas. 

The Pencil Code implements thermodynamic evolu¬ 
tion by solving for the entropy S = (InP - 7 lni 7 ). We 
use thermal relaxation for non-isothermal runs, driving 
the temperature to a reference temperature T^ef within a 
time T, according to 


T 


DS 


(T-T,ef) 


+ Ah/ 


( 8 ) 


where both Tj.ef ^ radially-varying. We take Aef 
to be the initial temperature and r proportional to the 
orbital period 


t = toHo/H, (9) 

where i? = is the Keplerian angular fre¬ 

quency. The subscript "0" refers to the values of the 
quantities at an arbitrary reference radius rg = 1. The 
reference timescale for thermal relaxation, Tq, is a free 
parameter of the model reflecting the radiative cooling 
timescale of the gas. The function 


nh = l^sh(V-M)" (10) 

is the heating due to dissipation of shocks. 

Sixth-order hyper-dissipation terms are added to the 
evolution equations to provide extra dissipation near 
the grid scale, as discussed in Lyra et al. (2008). These 
terms are needed for numerical stability because the 
high-order scheme of the Pencil Code formally lacks er¬ 
ror terms to produce stabilizing numerical dissipation 
(McNally et al. 2012). 

The simulations are run in the inertial frame, centered 
at the barycenter of the system. The positions of the star 


and planet are evolved on circular orbits by direct in¬ 
tegration using an N-body routine, employing the same 
3rd-order Runge-Kutta algorithm used for updating the 
gas grid. Since we are modeling disks with high-mass 
protoplanets, this will allow us to capture any dynam¬ 
ical effects on the gas induced by stellar wobble. The 
gravitational potential ^ acting on the star, planet, and 
gas is the sum of the potential contributions of the star 
and planet (i.e. the gravity of the gas is ignored): 


N 


E 

i 



( 11 ) 


where N = 2 , M/ is the mass of particle /, TZi = \r- r/1 is 
the distance between the zth particle and a gas parcel or 
massive particle for which the potential is being calcu¬ 
lated, and bf is the smoothing radius of the zth particle, 
which is used to prevent singularities. The gravitational 
smoothing radius of the planet is set to its Hill radius. 
The star does not require a smoothing radius because it 
lies outside the gas grid. 

For units, the planetary semi-major axis and orbital 
velocity are used as the units of length and velocity, re¬ 
spectively. This in turn implies 2n time units per orbit. 


2.2. Run parameters 

All simulations are performed on a cylindrical grid 
with {Nr, Nq) = (768,1024), with 0.4 < r < 12 (in units of 
the planet's semi-major axis). Grid points are spaced in 
r according to a power law such that Ar oc r^-^, providing 
higher resolution in the vicinity of the planet compared 
with the outer disk (for our choice of initial tempera¬ 
ture gradient, this non-equidistant grid corresponds to 
a constant number of radial grid cells per scale height). 
For all simulations, the initial density distribution is ax- 
isymmetric and decreases with the square root of the 
radius, with Nq = 1 (the physical units of surface den¬ 
sity can be chosen arbitrarily since the self-gravity of 
the gas has not been included). The sound speed at 
the reference radius was set to 5 x 10 “^, which corre¬ 
sponds to a temperature of -100 K, assuming a 1M© 
star with planetary semi-major axis 5.2 AU, 7 = 7/5, 
and mean molecular weight 2.4 (corresponding to a 5:2 
hydrogen/helium mixture by mass). The initial sound 
speed is set to 

Cs = (12) 

In order to study the role of thermodynamical evo¬ 
lution in disks containing high-mass protoplanets, we 
perform four sets of simulations. First, we perform sim¬ 
ulations for three planet-to-star mass ratios q = 10 “^, 
5 X 10“^, and 10“^, with locally isothermal and adia¬ 
batic disks. We run one adiabatic simulation with shock 
heating artificially turned off in the entropy equation 
(Ah = 0) / while keeping the acceleration due to shock 
viscosity turned on in the momentum equation. This 
is done to determine the role of shock heating in the 
non-barotropic simulations. For this test case, q = 10“^. 
Next, we run three simulations with different values of 
the reference thermal relaxation timescale Tq, also with 
q = 10“^. Finally, we run three additional simulations 
with mass ratio 10 “^, but using different radial power 
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TABLE 1 

Simulation parameters 


Run q 


7 

To 

Shock 

heating 


Mass, thermal evolution; Eigs. 1-4 

A 

10-3 

-1 

1 

0 

— 

B 

5 X 10“3 

-1 

1 

0 

— 

C 

10-2 

-1 

1 

0 

— 

D 

10-3 

-1 

1.4 

oo 

Yes 

E 

5 X 10-3 

-1 

1.4 

oo 

Yes 

E 

10-2 

-1 

1.4 

oo 

Yes 

G 

Viscous heating; Fig. 3 
10-2 -1 1.4 00 

No 


Relaxation time; Eig. 5 


H 

10-2 

-1 

1.4 

0.1 

Yes 

I 

10-2 

-1 

1.4 

1 

Yes 

J 

10-2 

-1 

1.4 

10 

Yes 

K 

Temperature power law; Eig. 6 

10-2 -0.5 1.4 00 Yes 

L 

10-2 

-0.2 

1.4 

00 

Yes 

M 

10-2 

0 

1.4 

oo 

Yes 


laws for the initial temperature (all adiabatic). Table 1 
contains run parameters for all simulations, thirteen in 
total, where q is the planet-to-star mass ratio, j 6 is the 
slope of the radial temperature power law, 7 is the adi¬ 
abatic index of the gas, and Tq is the cooling timescale 
in units of orbits at the reference radius. Locally isother¬ 
mal runs are denoted by Tq = 0 and 7 = 1 , while adia¬ 
batic runs are denoted by Tq = 00 . All simulations are 
run for 100 planetary orbits. 

3. RESULTS 
3.1. Global disk properties 
3.1.1. Planet mass and thermal evolution 

Fig. 1 shows the surface densities in the inner disk 
(r < 4) after 100 orbits for the locally isothermal (tq = 0) 
and adiabatic (tq = 00 ) runs and for three values of the 
planet-star mass ratio q. For all q, there is a clear differ¬ 
ence between the global morphologies of the isothermal 
and adiabatic disks. All three isothermal simulations 
show well-behaved spiral structure, a clear gap, and an 
accompanying vortex, as expected. In contrast, large- 
scale turbulence is clearly present in the q = 5 x 10 “^ and 
10 “^ adiabatic runs, in which the expected gaps and vor¬ 
tices appear to be completely missing. For the q = 10"^ 
adiabatic run, the presence of global turbulence is not 
as apparent, but the gap depth and the accompanying 
vortex are visibly diminished. 

Since the turbulence seen in Fig. 1 is present only in 
runs with non-barotropic equation of state, we hypoth¬ 
esize that it is the result of buoyant oscillations trig¬ 
gered by hot gas that has been tidally compressed by 
the planet. 

To test this hypothesis and confirm the originating re¬ 
gion of any such oscillations, we compare the results of 
each simulation to the Solberg-Hoiland stability crite¬ 
rion. The criterion for hydrodynamic stability is given 

by 

k^ + N^> 0, (13) 


where the epicyclic frequency 


K = 


■^_d 

r3 



(14) 


and the characteristic frequency of buoyant oscillations, 
the Brunt-Vaisala frequency. 


N = 


■ 1 dP , 

rids 

1 dpy 

T dr ' 

vT d7 

7 P dr / 


1/2 


(15) 


In the case of a steep entropy gradient, will become 
negative, leading to buoyant instability. For positive 
the velocity shear of the disk will tend to damp buoy¬ 
ant oscillations. We calculate A for each simulation at 
several timesteps, where 


A = l + 




(16) 


is the left-hand side of Equation (13), normalized by k^. 
For A < 0, the local shear of the disk will not stabilize 
buoyant oscillations. If the turbulence observed in Fig. 1 
is due to buoyant instabilities, then A will achieve nega¬ 
tive values only in the adiabatic simulations. This is be¬ 
cause in the locally isothermal simulations neither the 
temperature profile (which remains constant through 
the simulation) nor the density fluctuations are suffi¬ 
cient to give rise to a steep entropy gradient. 

For the locally isothermal simulations, we indeed find 
that A never becomes negative. In the adiabatic case, on 
the other hand, for q = 10“^ (run F), A, shown in left 
panel of Fig. 2, reaches negative values after less than 
one orbit in the region of the outer spiral wake of the 
planet, as well as, to a lesser extent, the inner spiral 
wake. After 100 orbits, zones of buoyant instability are 
seen throughout the disk, as shown in the right panel of 
Fig. 2. For the q = 10“^ case (run D), A achieves nega¬ 
tive values for the first ~15 orbits, then stays positive for 
the remainder of the simulation, explaining the relative 
similarity in resulting surface densities between the cor¬ 
responding isothermal and adiabatic cases (runs A and 
D, respectively). These results confirm that buoyant in¬ 
stabilities are consistent with the large-scale disruption 
of the usual planet-disk interaction features seen in the 
lower panels of Fig. 1. 


3.1.2. Shock heating 

In order to determine the role of shock heating in gen¬ 
erating the observed buoyant instability, we run another 
simulation with an adiabatic disk and a 10 Mj planet 
where we exclude heating of the gas due to shocks 
(Tgh = 0; run G). The resulting surface density after 100 
orbits, shown in the right panel of Fig. 3, resembles that 
of the locally isothermal case, where the gap and vor¬ 
tex are visible. This test conclusively shows that the en¬ 
tropy generated in shocks is responsible for the novel 
phenomena described in the present work. 


3.1.3. How strong are the shocks? 

In Fig. 4 we show radial profiles of temperature (nor¬ 
malized by the initial temperature profile) and Mach 
number Ma for 100 orbits in the adiabatic, 10 Mj case 
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Fig. 1.— Surface densities in the inner disk (r < 4) after 100 orbits for locally isothermal (tq == 0) runs A, B, and C, and adiabatic (tq = oo) runs D, 
E, and R Turbulence due to buoyancy appears in adiabatic runs and increases with q (masses based on M star ~ 1M©). 



Fig. 2.— A for q = 10“^ after 1 orbit (left panel) and 100 orbits (right 
panel) for the inner disk (r < 4) in the adiabatic case (run F). Buoyant 
instability occurs for A < 0. It begins in the outer spiral wake pro¬ 
duced by the planet after 1 orbit, while instability is triggered in the 
inner spiral wake and outer disk soon after. 

(run F) at the azimuthal position of the planet. Shocks 
propagate throughout the disk for the duration of the 
simulation; however, they do not grow rapidly and 
monotonically like the temperature, with the strongest 
shocks (Ma > 1.5) subsiding after ~ 50 orbits. The tem¬ 
perature increases steadily throughout the disk over the 
course of the simulation. This shows that the tempera¬ 
ture increase is not due to a single shock, but due to the 


entropy injection of multiple shocks, over the course of 
many orbits. 

3.1.4. Thermal relaxation timescale 

To investigate the importance of cooling in non- 
barotropic disks with high-mass planets, we run several 
simulations with q = 10"^ for cooling times Tq = 0.1, 1, 
10, and 00 orbits (runs H, 1, J, and F, respectively). Sur¬ 
face densities after 100 planetary orbits are shown in 
Fig. 5. Not surprisingly, simulations with shorter relax¬ 
ation times (0.1 and 1 orbits) more closely resemble the 
isothermal case: buoyant instabilities and subsequent 
turbulence are minimal, and so the gap and vortex are 
preserved. For longer cooling timescales (10 orbits and 
adiabatic), the shock heating drives buoyant instabil¬ 
ities, generating turbulence that dominates the global 
structure of the disk, disrupting the gap and vortex. 

Because the wake is stationary in the reference frame 
of the planet, a parcel of gas will be excited by one of the 
spirals once per synodic period Tsyn (with respect to the 
planet). Therefore, for Tq < Tsyn, the parcel will cool ra- 
diatively before encountering the next shock, and thus 
the energy contributions of successive shocks will not 
continue to increase the temperature of the parcel be¬ 
yond that of a single shock. For Tq > Tsyn, however, the 
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Fig. 3.— Adiabatic simulation of planet-disk interaction of a lOMj planet (left) with shock heating and (right) without shock heating. This test 
conclusively shows that the energy dissipated in shocks is driving the instability. 



Fig. 4.— Temperature and Mach number for the inner disk (r < 5) of 
run F, at the azimuthal position of the planet. Heating due to shocks 
initially occurs close to the planet, but spreads to the rest of the disk 
over the course of the simulation. 

gas parcel will not have time to cool to its original tem¬ 
perature before encountering the next shock, therefore 
the temperature of the gas will continue to rise, even 
though the energy contributions of each shock will be 
the same as in the low Tq case (since Eq. 10 does not de¬ 
pend on cooling timescale). The emergence of the ob¬ 
served buoyant instabilities therefore depends on the 
relative orbital and cooling timescale. Its strength de¬ 
pends additionally on the energy contribution of each 
shock. These are a function of the characteristic differ¬ 
ential velocity of the gas flows associated with the spi¬ 
ral wakes, which is in turn determined by the planet's 
mass. 

3.1.5. Temperature power law 

Because the buoyant instability results from the pres¬ 
ence of a strong entropy gradient, we compare sim¬ 
ulations with differing initial temperature power law 
slopes j6 to determine whether a steeper initial entropy 
gradient enhances the buoyant instabilities generated 


To = 0.1 orbits tq = 1 orbits 



X 


Fig. 5.— Surface densities for the inner half of the disk (r < 6) after 
100 orbits for four different cooling timescales (runs H, I, J, and F), for 
a planet with q = 10“^. 

locally by the planet. We once again use q = 10“^ and 
an adiabatic disk, and conduct runs with j6 = 0, -0.2, - 
0.5, and -1, corresponding to runs M, L, K, and F, re¬ 
spectively. The value for run L was chosen to flatten 
the initial entropy gradient, i.e., ds/dr = 0. For = 0, 
A is initially positive throughout the disk, while it be¬ 
comes more negative for decreasing j6, so that for these 
values the disk is initially slightly Solberg-Hoiland un¬ 
stable (see Eq. 16). 

The resulting surface densities after 100 planetary 
orbits for these simulations are shown in Fig. 6. A 
steep temperature slope (e.g., j6 = -1) appears to lead 
to stronger turbulence, though significant turbulence is 
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Fig. 6 .— Surface densities for the inner half of the disk (r < 6) after 
100 orbits for four different temperature power law slopes (runs M, L, 
K, and F). Run M is initially globally stable according to the Solberg- 
Hoiland criterion, but still experiences instability due to heating from 
the planetary wakes (see Fig. 2). 

still observed for j6 = 0. The enhancement of the tur¬ 
bulence when the disk is globally unstable supports our 
hypothesis that the planet generates turbulence by driv¬ 
ing local regions into instability, even when the disk is 
globally stable. 


3.2. Angular momentum transport 

In order to examine angular momentum transport in 
the disks simulated in this work and thereby identify 
possible implications for disk evolution, we calculate ef¬ 
fective Shakura & Sunyaev (1973) A:-viscosities over time 
for simulations of locally isothermal and adiabatic disks 
with q = 10“^ (runs C and F). We define an effective vis¬ 
cosity V in terms of oc: 

v = 0 Lc\in, (17) 

which drives an accretion rate (assuming a steady state) 


m = -27TvE 


din f2 
dlnr ' 


= 3/11/27. 


(18) 


We define also the oc value as a function of the kinetic 
stress SUfSucp (where 5ui = Ui - (zi/), and (X) represents 
the azimuthal average of X) 


2 {SUy SU(p^ 

( 4 ) 


(UrU^) - {Ur){u^) 


(19) 
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Fig. 7.— Azimuthally-averaged effective a-viscosities for locally 
isothermal (C; left) and adiabatic (F; right) runs, over the course of 
100 planetary orbits, for r <9. The high oc values in the isothermal 
case, taken at face value, imply nearly supersonic turbulence, and in¬ 
dicate that the local a-disk prescription does not properly describe the 
disk's evolution and that the thin-disk approximation may be break¬ 
ing down. The thin-disk approximation also limits the adiabatic case: 
in three dimensions the shocks will be weaker as vertical expansion is 
allowed. The results nevertheless highlight that shocks in the wakes 
of high-mass planets are important for mass redistribution. 

case (left panel of Fig. 7) the oc values are large (reach¬ 
ing 0.5) but occur close to the planet and gap, while 
the disk beyond r =2-3 remains relatively quiescent. By 
contrast, the turbulence in the adiabatic case propagates 
outward through the disk throughout the entire simula¬ 
tion. The oc in this case are lower, reaching 0.05, an order 
of magnitude less. The globally-averaged oc between 60- 
100 orbits is ^0.01-0.02. The relatively large value of oc 
indicates rapid inward motion of the gas. Defining the 
accretion timescale 

Tacc= m/rh ~ f2r^{cloc)~^, (20) 

we find that for a lOMjup planet orbiting a 1M© star at 
5.2 AU, the characteristic accretion timescale of the disk 
is on the order of 10^ years. 

We caution that this particular result pertaining to ac¬ 
cretion is a function of the cooling time, as for grav¬ 
itational instability (Durisen et al. 2007; Meru & Bate 
2011, 2012), or any shock-driven instability. The 
high oc values in the isothermal case (~ 0.5) if taken at 
face value, imply that the rms velocity of the turbu¬ 
lence is approaching supersonic values (z/rms ~ Cg ~ 
0.7cs). However, this is mainly an indication that the 
local A:-disk prescription does not properly describe 
the disk's evolution and that the thin-disk approxima¬ 
tion may be breaking down. In the adiabatic case, 
even though the effective oc values are reasonably sub¬ 
sonic, the large-scale turbulence is similarly implicated. 
In three-dimensional runs we expect the shocks to be 
weaker because of the extra degree of freedom (verti¬ 
cal expansion such as in shock bores; Boley & Durisen 
2006) and the efficient radiative cooling of the disk's up¬ 
per layers, even though the midplane is optically thick 
and essentially adiabatic. 

3.3. BoxzyHydro simulations 


Fig. 7 shows oc for the aforementioned simulations, 
calculated once per orbit for 100 orbits. 

We find that the evolution of the disk driven by the 
effective oc viscosity differs substantially between the lo¬ 
cally isothermal and adiabatic cases. For the isothermal 


Muller & Kley (2013) perform 2D simulations of 
disks containing protoplanets up to 16 MJup using 
the FARGO code (Masset 2000; Baruteau 2008), using 
isothermal and radiative viscous accretion disk models. 
Their results do not show the strong temperature in- 
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Fig. 8 .— BoxzyHydro simulations with a 10Mj planet. Surface density, temperature, and pressure are all normalized to the initial values at 
the planet's semi-major axis. In the upper panels we show an isothermal control run. A clear spiral structure is seen, along with edge waves by 
the outer gap wells, evidencing the initial stages of vortex formation. In the lower panels we show an adiabatic run, in which these features are 
gone. The shock-driven buoyant instability is reproduced with an independent code. BoxzyHydro solves for the equations of hydrodynamics in 
conservative form and uses an approximate Riemann solver for determining fluxes at cell interfaces. It does not require artificial viscosity terms in 
the energy equation for shock capturing. This rules out the possibility that the instability is effected by Pencil's numerical shock capturing scheme, 
evidencing a physical origin for the adiabatic case. 


crease seen here. We first suspected that the difference 
lies in the applied cooling laws: while Muller & Kley 
(2013) use a cooling function that depends on p and T, 
we apply a fixed, radially dependent cooling. To check 
this, we performed a calculation defining a dynamical 
cooling time depending on p and T as well. The results 
were consistent with the fixed, radially dependent cool¬ 
ing time. Unfortunately, Muller & Kley (2013) do not 
specify whether their code included artificial viscosity 
in order to properly treat shock heating. If they did, 
it would enter in their viscous energy dissipation term 
Qvisc/ but its form was not specified. A meaningful com¬ 
parison between our work and theirs is therefore unfea¬ 
sible. 

To resolve this disagreement, and confirm the occur¬ 
rence of buoyancy-induced turbulence in a 2D, adia¬ 
batic disk with a massive protoplanet, we reproduce run 
F (10 Mj, adiabatic) using an independent code, Boxzy¬ 
Hydro (Boley et al. 2013). BoxzyHydro solves the 
equations of hydrodynamics in conservative form on a 
Cartesian grid. An approximate Riemann solver is used, 
along with standard flux limiting, to describe the fluid 
states at cell interfaces (e.g., Toro 2009). This allows the 
code to capture shocks without including explicit artifi¬ 
cial viscosity terms in the gas equations. 

Using BoxzyHydro, we model the disk within the 
radial range 0.4 < r < 4 (where the planet's semi-major 
axis ro = 1, as before). The initial setup uses the same 


density and temperature profiles, differing only in the 
details of the boundary conditions and the smoothing of 
the planet gravitational potential. As the code is Carte¬ 
sian, there is no strict inner boundary for the grid. The 
region interior to r = 0.3 (~1.6 AU for rg = 5.2 AU) 
is locally isothermal to avoid intense heating close to 
the star. The planet potential is piece-wise smoothed 
with a cubic spline inside the Hill radius, while keeping 
the Newtonian profile outside. In comparison, in the 
Pencil simulations, the Plummer smoothing used there 
(Eq. (11)) results in a shallower potential well immedi¬ 
ately outward of the Hill radius, therefore slightly un¬ 
derestimating the strength of the shocks. 

We run the BOXZYHYDRO simulation for 21 orbits. 
The result is shown in Fig. 8. The upper panels show an 
isothermal control run, while the lower panels show the 
adiabatic simulations. The isothermal run shows well- 
behaved spiral structure, and edge waves in the outer 
region of the gap, produced by the growing phase of 
the RWI. In the adiabatic run these features are gone, 
and instead the code reproduces the main features seen 
in Pencil's run F at the same time. At this point, turbu¬ 
lence has already begun to proliferate through the disk, 
diminishing the gap (which in run F disappears com¬ 
pletely at f ~ 50 orbits), and suppressing the formation 
of both inner and outer vortices. These results rule out 
the possibility that the observed turbulence is an artifact 
of the shock capturing scheme used in the Pencil Code. 
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The confirmation by an independent code strengthens 
the case that the result is not a numerical artifact. 

4. DISCUSSION AND CONCLUSIONS 

Using 2D, global hydrodynamics simulations of disks 
with embedded massive planets, we have shown that 
shocks generated by high-mass planets can significantly 
raise the temperature of gas in disks. While any given 
shock is typically weak, with Mach numbers of order 
unity for 5Mj and up to 4-5 for 10 Mj, the temperature 
increase is due to the cumulative effects of shocks, as or¬ 
biting gas repeatedly meets the spiral wake (which re¬ 
mains stationary in the reference frame of the planet). 
If the gas can cool within a synodic period, the energy 
radiates away between shocks and the situation resem¬ 
bles the isothermal case. If, on the other hand, the cool¬ 
ing time exceeds the synodic period, the mean tempera¬ 
ture increases in time. We then see buoyant instabilities 
driven by both the inner and outer spiral wakes, cre¬ 
ating sustained, large-scale turbulence that significantly 
alters the global structure of the disk. 

These results are of particular significance to the fre¬ 
quent interpretation that spiral morphologies are sign¬ 
posts of embedded protoplanets in transition disks. Rts 
to observed spiral structures are often attempted (e.g., 
Muto et al. 2012) by making use of linear spiral den¬ 
sity wave theory, which depends on the temperature 
gradient in the disk, as well as the actual disk tempera¬ 
tures (disk aspect ratio). However, because this theory 
is linear, it relies on assumptions of low planetary mass 
and local isothermality. In the past few years, some spi¬ 
ral features that have been observed that present prob¬ 
lems in interpreting them as having a planetary origin, 
if linear theory is employed. For example, Juhasz et al. 
(2014) found that the observed spirals in scattered light 
would require an increase by a factor of either > 3.5 in 
surface density or ^0.2 in pressure scale height above 
the background disk to generate them. They favor the 
increase in pressure scale height to reproduce the ob¬ 
servations. This is consistent with the observed spirals 
showing in general an opening angle wider than ex¬ 
pected from the background temperature. Spiral fitting 
to the disk MWC 758 (Benisty et al. 2015) at radii r p^^SO- 
150 AU needs a large aspect ratio, corresponding to to a 
disk temperature of 300 K at 55 AU, while the surround¬ 
ing gas is only at p:i50K. Finally, recent observations of 
HD 100546 (Currie et al. 2014) show a planetary can¬ 
didate without an associated spiral, and a feature that 
resembles a spiral arm, but at 90° away from the can¬ 
didate planet. Moreover, the spiral feature shows lit¬ 
tle polarization, implying thermal emission at roughly 
1000 K. These features are a challenge to linear spiral 
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density theory but can easily be explained in light of 
our model, which predicts a significantly different qual¬ 
itative behavior for spiral structures generated by high- 
mass planets in radiatively inefficient disks. 

In the following paper in this series, we will describe 
three-dimensional simulations to further explore the 
role of shock heating from the wakes of massive planets 
in the structure and evolution of the disk. We caution 
again that in three dimensions the results will likely dif¬ 
fer. Shocks should be weaker for the same planetary 
mass, given the possibility of vertical expansion. It is 
possible that the turbulence seen in this work would be 
seen in another region of parameter space; for instance, 
it may require a higher protoplanet mass to achieve an 
effect comparable to that observed in the current work. 
Moreover, even though the midplane of the disk is opti¬ 
cally thick and close to adiabatic, the upper layers might 
cool efficiently. In this case, the diffusion timescale in 
the vertical direction may control a given disk's evolu¬ 
tion. We note, however, that the buoyant instability may 
share similarities to other shock-driven processes, such 
as gravitational instability (Durisen et al. 2007), where 
reasonable agreement has been shown between 2D and 
three-dimensional simulations (see Rice et al. 2014 and 
references therein). Evidence for similar qualitative be¬ 
havior is seen in Boley & Durisen (2006), who modeled 
shocks due to 2.5 Mj planets, observing shock bores as 
gas accelerates upwards and breaking waves as the gas 
descends back onto the disk, again generating turbu¬ 
lence, though by a different mechanism. 

We finally note that heating and buoyancy-induced 
turbulence may also affect planet migration. Investigat¬ 
ing this effect (as in Zhu et al. 2012, but for high-mass 
protoplanets) would require simulations to be carried 
out over many more orbits, and for the gravitational 
back-reaction of the gas on the planet to be included. 
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